%{
    This file is part of StemCellQC, a video bioinformatics software
    toolkit for analysis of phase contrast microscopy videos.
    Copyright 2013-2015 Vincent On. [vincenton001-at-gmail.com]

    StemCellQC is free software: you can redistribute it and/or 
    modify it under the terms of the GNU General Public License as 
    published by the Free Software Foundation, either version 3 of the 
    License, or (at your option) any later version.

    StemCellQC is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with StemCellQC.  If not, see <http://www.gnu.org/licenses/>.
%}

function skip_frame = fun_automate( index , seg_method, segmentation_gain, skip_frame, fudgeFactor )
global mov

try
    switch seg_method
        case 1
            filtSize = 1;
            
            % Normalize the image
            temp = double(mov(index).data) - min(min( double(mov(index).data) ));
            temp = temp ./ max(max( temp ));
            
            % Blur it a little bit
            temp = imfilter( temp, ...
                fspecial( 'gaussian', (1+2*2*filtSize) ), ...
                'symmetric' );
            
            % Entropy filter: The cell is a uniform texture
            temp = entropyfilt( temp, true((1+2*filtSize)) );
            temp = temp - min(min( temp ));
            temp = temp ./ max(max( temp ));
            % I_min = min(min( double(mov(index).data) ));
            % I_max = max(max( double(mov(index).data) ));
            % temp = ((double(mov(index).data) - I_min)./(I_max - I_min))*255;
            
            % Opening: we want small 'balls' touching the cell to stop touching
            %{
        TODO: The following section could be improved. It doesn't always
        prevent non-cell objects from being counted as a part of the cell.
            %}
            temp = imopen( temp, ...
                strel('disk',round((1+2*3*filtSize)*segmentation_gain)) ); %try larger disc sizes if segmentation is picking up too much
            
            % Now threshold
            [ thres,EM] = graythresh( temp );
            temp = im2bw( temp, thres );
            
            % temp = imopen( temp, ...
            %     strel('disk',(1+2*3*filtSize)) );
            
            % Retain only the largest object
            [Largest, ~] = fun_returnLargestObj( temp );
            
            % Check if the largest object is in the correct region
            if skip_frame == 0
                if index > 1;
                    stop = 0;
                    iter = 1;
                    
                    while stop == 0;
                        
                        centroid = regionprops(Largest, 'Centroid');
                        cx = centroid.Centroid(1);
                        cy = centroid.Centroid(2);
                        
                        [t,b,l,r] = bound_box(mov(index-1).features{3},mov(index-1).features{4},150);
                        
                        if cx < l || cx > r || cy < t || cy > b
                            temp(Largest~=0) = 0;
                            [Largest, ~] = fun_returnLargestObj( temp );
                            
                        else
                            stop = 1;
                            
                        end
                        iter = iter + 1;
                        
                        if iter > 10;
                            stop = 1;
                        end
                    end
                    
                    %         %Check area of previous segmentation
                    area = regionprops(Largest, 'Area');
                    prev_area = mov( index - 1 ).features{ 1 };
                    if area.Area < ( .5 * prev_area ) || area.Area > ( 2 * prev_area )
                        error('Change in Area is too large');
                    end
                end
            end
            
        case 2
            I = mov(index).data;
            
            % Detect Cell Edges
            [~, threshold] = edge(I, 'sobel');
%             fudgeFactor = .5;
            BWs = edge(I,'sobel', threshold * fudgeFactor);
            
            % Dilate the edges
            se90 = strel('line', 3, 90);
            se0 = strel('line', 3, 0);
            BWsdil = imdilate(BWs, [se90 se0]);
            
            % Fill
            BWdfill = imfill(BWsdil, 'holes');
            
            % Smooth Objects
            seD = strel('diamond',1);
            BWfinal = imerode(BWdfill,seD);
            BWfinal = imerode(BWfinal,seD);
            temp = BWfinal;
            [Largest, ~] = fun_returnLargestObj( temp );
    end
    temp = Largest;
    % Fill it
    temp = imfill( temp, 'holes' );
    
    % SAVE IT HERE
    mov(index).segmented = temp;
    clear temp
    
    %====
    %Part 2: gettting region statistics
    
    stats = regionprops( mov(index).segmented, ....
        'Area', ...
        'Perimeter', ...
        'Centroid', ...
        'Extent', ...
        'Solidity', ...
        'Orientation', ...
        'MajorAxisLength', ...
        'MinorAxisLength', ...
        'Eccentricity', ...
        'Extrema' );
    
    % compute statistics perimeter pixels with respect to the centroid
    D = cen2peri_dist( mov(index).segmented );
    Min_rad = min(D);
    Max_rad = max(D);
    Avg_rad = mean(D);
    
    % compute difference of average intensities of background and silhouette
    % regions
    [avg_intensity, max_intensity, min_intensity, bright_area] = intensity_extract(index, mov(index).segmented, 75);
    
    bright_ratio = bright_area/stats.Area;
    
    % compute number of protrusions and area of protrusions
    [N_protrusions, Prot_area] = protrusions(index,mov(index).segmented,25);
    Prot_ratio = Prot_area/stats.Area;
    
    mov(index).features( 1 ) = { stats.Area };
    mov(index).features( 2 ) = { stats.Perimeter };
    mov(index).features( 3 ) = { stats.Centroid(1) };
    mov(index).features( 4 ) = { stats.Centroid(2) };
    mov(index).features( 5 ) = { stats.Extent };
    mov(index).features( 6 ) = { stats.Solidity };
    mov(index).features( 7 ) = { stats.Orientation };
    mov(index).features( 8 ) = { stats.MajorAxisLength };
    mov(index).features( 9 ) = { stats.MinorAxisLength };
    mov(index).features( 10 ) = { stats.Eccentricity };
    mov(index).features( 11 ) = { stats.Extrema };
    mov(index).features( 12 ) = { Min_rad };
    mov(index).features( 13 ) = { Max_rad };
    mov(index).features( 14 ) = { Avg_rad };
    mov(index).features( 15 ) = { avg_intensity };
    mov(index).features( 16 ) = { max_intensity };
    mov(index).features( 17 ) = { min_intensity };
    mov(index).features( 18 ) = { bright_ratio };
    mov(index).features( 19 ) = { N_protrusions };
    mov(index).features( 20 ) = { Prot_ratio };
    
    skip_frame = 0;
catch
    skip_frame = skip_frame + 1;
    disp( ['Skiping frame ' num2str( index )]);
end